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Slow relcixation in the two dimensional electron plasma under the strong magnetic 
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We study slow relaxation processes in the point vortex model for the two-dimensional pure 
electron plasma under the strong magnetic field. By numerical simulations, it is shown that, 
from an initial state, the system undergoes the fast relaxation to a quasi-stationary state, and 
then goes through the slow relaxation to reach a final state. From analysis of simulation data, 
we find (i) the time scale of the slow relaxation increases linearly to the number of electrons 
if it is measured by the unit of the bulk rotation time, (ii) during the slow relaxation process, 
each electron undergoes an superdiffusive motion, and (iii) the superdiffusive motion can be 
regarded as the Levy flight, whose step size distribution is of the power law. The time scale 
that each electron diffuses over the system size turns out to be much shorter than that of the 
slow relaxation, which suggests that the correlation among the superdiffusive trajectories is 
important in the slow relaxation process. 

KEYWORDS: non-neutral plasma, two-dimensional turbulence, numerical simulation, long range force, 
non-extensive system, slow relaxation, anomalous diffusion. Levy flight 
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1. Introduction 

The system of pure electron plasma under the strong 
magnetic field has been of scientific interest for various 
aspects; Such a system does not satisfy the basic premise 
of the Boltzmann statistics, i.e. the existence of a subsys- 
tem weakly coupled to the rest, because the interaction 
between the constituent particles is long-range, conse- 
quently any part of the system interacts strongly with the 
rest. Under a certain experimentally realizable condition, 
the system behaves like a two-dimensional (2-d) system 
and is shown to be described as a 2-d point vortex sys- 
tem with a unit vorticity, whose continuum description 
reduces to the 2-d Eulcr equation of the incompressible 
inviscid fluid. ^ 

Historically, the point vortex system is the system for 
which Onsager developed the idea of negative tempera- 
ture in his statistical theory and predicted the equilib- 
rium states of the large vortex cluster in the relaxation of 
2-d high Reynolds number fluid. His theory has been 
advanced, ''^^ and has been examined by simulations of 
several hundred particles (or vortices) to show that the 
equilibrium state is the vortex cluster that maximizes the 
one-body entropy.^ 

In contrast with these results, in recent 2-d electron 
plasma experiments, the system has been demonstrated 
to relax into several kinds of "stationary" or "metaequi- 
librium" states, including the minimum enstrophy state^ 
and the vortex crystal state. ^ These are quite differ- 
ent from the above equilibrium state. There have been 
applied several statistical theories to understand these 
relaxed states, which includes the minimum enstrophy 
theory,^ the fluid entropy theory, the Tsallis statis- 
tics,^'^ the regional entropy theory for the vortex crys- 
tal,^'* etc., but none of these theories have successfully 
described these states on a general ground. 
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These two contradicting results suggest that the sys- 
tem is trapped in some states before reaching the equi- 
librium in the experiments. In fact, we show in the 
present work that these stationary states are only quasi- 
stationary and relaxes into the equilibrium state very 
slowly. This process of slow relaxation is the major sub- 
ject of the present paper. 

Actually, for the stellar system with gravitation, which 
is another example with long-range interaction, it has 
been known that there are two types of relaxation: the 
violent relaxation and the coUisional relaxation. The vi- 
olent relaxation^^' is the fast process which is caused 
by the complex mixing due to the macroscopic collec- 
tive interaction. The dynamics is described by the Vlasov 
equation and the two-body "collisions" do not play any 
role. In contrast, the coUisional relaxation^°'^^ is the slow 
process caused by the collisions of the individual particles 
and the relaxation time scale is much longer than that 
of the violent relaxation. Such a relaxation has already 
been found and analyzed in the molecular dynamics sim- 
ulations of self-gravitating system. The separation of 
these two relaxations is considered to be peculiar to the 
systems with long-range interaction. 

In the present system of the 2-d electron plasma under 
the strong magnetic held, there are also a fast and a slow 
relaxation. The fast relaxation from an initial state takes 
place through the violent motion of the stretching and 
folding in the density field due to the macroscopic collec- 
tive motion, and leads the system to a quasi-stationary 
state. The dynamics is described by the 2-d Euler equa- 
tion and the individual two-body collisions do not play 
any role. The typical time scale for this fast relaxation is 
the rotation time of the bulk of the system Trot . 

In contrast, the relaxation after the quasi-stationary 
states towards the equilibrium is much slower than the 
fast relaxation. 

Chavanis has presented a theoretical analysis on the 
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dynamics of this system^'*' using the analogous idea 
developed in the sclf-gravitating system. He predicted 
that the ratio of the slow relaxation time to the fast 
relaxation scales with the total number of particles N as 
N/ InN, based on the idea that the relaxation takes place 
through the diffusion of the individual particles caused by 
collisions, which can be ignored in the violent relaxation. 

For the neutral plasma under the strong magnetic field, 
the theoretical analysis by Taylor and McNamara^^ has 
shown that the individual particle motion converges to 
the normal diffusive behavior only slowly, and the lim- 
iting value of the diffusion coefficient have been found 
to depend on the system size, which implies anomalous 
diffusion is taking place. Dawson et al. have pointed out 
that such behavior of diffusion is attributed to occasional 
long jumps of the particles convected by large vortices. 

In some experiments of the 2-d fluid with vortices, the 
diffusion of tracer particles have been often found to be 
anomalous, ^^'■^^ and a Levy flight analysis has been pro- 
posed.'^" 

In this paper, we present the results of the numerical 
simulations on the point vortex system. We observe the 
iV-dependencc of the fast and slow relaxation times and 
find that they are consistent with those estimated by 
Chavanis. Special attention is on the diffusion process 
during the slow relaxation, for which we find an anoma- 
lous diffusion and analyze it in terms of the Levy flight. 

The paper is organized as follows: In §2, we introduce 
the 2-d electron plasma system under a magnetic field 
and define the model. In §3, we briefly review some sta- 
tistical theories. The simulation method is explained in 
§4 and results are described in comparison with the sta- 
tistical theories in §5. We analyze these results in terms 
of the Levy flight in §6, and the summary and discussion 
are given in §7. 

2. Model system and its behavior 

2.1 Equations of motion 

The physical system we consider is the pure electron 
plasma in a cylindrical container with the strong mag- 
netic field B applied along the axis (Fig. 1). In a cer- 
tain situation,^' ^ it has been shown that the interaction 
between the electrons can be approximated by the two- 
dimensional Coulomb force and the drift velocity Vi of 
the i-th electron in the plane perpendicular to B is given 

by 



E(r,) X B 



(1) 



disregarding the cyclotron motion. Here E{r) is the elec- 
tric field and is the position of i-th electron in the 
plane. Thus the electrons do not repel each other and 
can be contained in the cylinder (Malmberg Trap). 

After some normalization, it can be shown^^ that the 
equations of motion for the electron position = {xi,yi) 
are given by 

dxi dH dyi dH 
dt dyi ' dt dxi ' 
with the Hamiltonian H , 



(2) 



N N N 



1 ^ 

iE 



(3) 



where G{ri^ Tj) is the 2-d Green function for the electric 
potential that satisfies V^G(r,r') — 5{r — r') with an 
appropriate boundary condition, Gm{r,r') is the elec- 
tric potential at r brought about by the mirror charge 
induced by the charge at r', and (j)i{r) is the electric po- 
tential due to all the electrons except for the ith one. 
In the present case, we consider the system which is in a 
cylindrical container with the radius Rw , then 0,; is given 

by31,32 
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(4) 



where z = x + and z* is the complex conjugate of 
z; the second term corresponds to the potential by the 
mirror charges. 

2.2 Euler equation 

Note that the density field 



i(r) = ^<5(r-r,(0) 



(5) 



with ri{t) being a solution of eq. (2) gives a singular 
solution for the set of equations, 

Dn _ dn 
iDt " 'dt 



v-Vn^ 0, 



(6) 



d(j> d(j) 
dy ' dx 



(7) 




0(V) 



Fig. 1. Schematic diagram of Malmberg trap of the pure electron 
plasma. 



= n, (8) 

where z denotes the unit vector perpendicular to the 
plane and V is the 2-d nabla. 

It can be shown that the density field n is equal to the 
vorticity w(r) = (V x v)z of the 2-d velocity field v and 
the velocity field v is solenoidal ( V • ti = 0) , therefore, the 
set of equations (6) — (8) is equivalent to that of the Euler 
equation for the 2-d incompressible inviscid fluid with a 
free-slip (no-stress) boundary condition, but the vorticity 
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i(r) = 



takes only a positive value in the present system. 

To see the basic motion and the time scale, we consider 
the "pancake" state with the density field, 

na (r < R) 
(r>i?) ' 

which is a steady solution of eq. (6) — (8). The system 
shows a rigid rotation with the period 

r An , ^ 

Trot = - — . (10) 

Vi-ot[r) no 

This period gives the approximate time scale for the 
macroscopic dynamics of the system of point vortices 
with the average density no « nN/B? . 

If the system contains macroscopic numbers of elec- 
trons that follow eq. (2), and their distribution can be 
represented by a smooth density field ny^{x,y) obtained 
through averaging over a finite size mesh, i.e. coarse 
graining, then the density field should follow the follow- 
ing partial differential equations: 



tem, 



Dnm 
Dt 



dr 



dt 



X V0n 



C(nm), 



dy ' dx 



(11) 



(12) 



= n„,. (13) 

These are almost the same with eqs. (6) — (8), but there 
is the extra term C{n^) in the right hand side of eq. 
(11). This comes from the averaging on n and represents 
microscopic processes, thus it is called the collision term, 
since it is analogous to the corresponding term in the 
Boltzmann equation for the kinetic theory of gases. Note 
that, in experimental situations, we observe the parti- 
cles density in a finite resolution, therefore the averaging 
operations are essential in macroscopic observations. 

When the system contains macroscopic number of par- 
ticles, the collision term is small, but it eventually causes 
the slow relaxation that leads the system to the equilib- 
rium state, and estimating the effects of the collision term 
is still under intense discussion. 



N 



(frr'^n{r) 
<fr (r X v{r))^, 



(15) 



is also a constant of the dynamics. 

Since the 2-d Euler equation (6) — (8) represents an 
area-preserving dynamics, it has a set of conserved quan- 
tities, called Casimir constants. 



Z= / <frf{n{r)), 



(16) 



where /(n) is an arbitrary function of the density field 
n. In the point vortex system, however, the Casimirs Z 
cannot be defined except for the linear combination of 
the total number of particles N 



(fr n{r) 



(17) 



since the density field n has singularity of the delta func- 
tion as in eq. (5). 

2.J^ Stability and evolution of states 

In this subsection, we focus on the macroscopic be- 
havior of the particle density n{r) obtained from coarse- 
grained observation (From now on, we drop the subscript 
'm' from the coarse-grained field, n,„(r) etc.). 

It can be shown that the rotationally symmetric state 
with the density n(r) being a decreasing function of r 
is not linearly unstable and is a numerically stable solu- 
tion of eqs. (6) — (8). By contrast, the ring state, where 
the electrons are distributed in a ring-shaped region, is 
linearly unstable, and its instability is called diocotron 
instability."^^ 

After the diocotron instability, nonlinear effects in- 
cluding the complex stretching and folding of the den- 
sity structure take places. These deformation causes a 
fine filament-like structure, which will soon become too 



2.3 Constants of the dynamics 

The Hamiltonian i? is a constant of the dynamics, and 
is expressed in terms of the field variables as 

H = f d^r f d^r' n{r)n{r')G{r,r') 



d^r n{r)(j){r) 



jd^r\v\rl 



which corresponds to the total energy. 

In the case of a system with circular symmetry, the 
total angular momentum / around the center of the sys- 
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Fig. 2. Time evolution of electron density distributions from (a) 
a single ring initial state, and (b) a double ring initial state. Ini- 
tial distributions are shown in the left-most figures and quasi 
stationary states are shown in the right-most figures. The con- 
ductor wall boundary is expressed by the solid arcs at the corner 
of the plotted area. 
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small to be observed in a finite resolution, and the system 
quickly approaches a quasi-stationary state. 

Figure 2 shows two examples of the time sequences, 
which start from the two unstable initial states and lead 
to two types of quasi-stationary states. Figure 2(a) shows 
the sequence that starts from a single-ring configuration, 
destabilizes into the mode three, and eventually falls into 
a singly peaked stable density distribution. In the case 
of Fig. 2(b), the instability of a double-ring initial con- 
figuration leads to a higher mode to break into many 
vortices, which merge occasionally while they undergo 
a collective chaotic motion, and eventually several sur- 
viving vortices form a regular structure, which docs not 
undergo further change during our simulation time; such 
a state is called a vortex crystal state. ^ 

These quasi-stationary states are in general not the 
equilibrium states;^^ the quasi-stationary states strongly 
depend on the initial states, which indicates that the 
system lacks the ergodicity in this time scale. For longer 
time scale, these quasi-stationary states relax into the 
equilibrium states, which we study in detail in the fol- 
lowings. 

3. Statistical theories for diffusion and relax- 
ation 

In this section, we briefly review some of the elements 
of statistical theory for the slow relaxation. 

3.1 Maximum one-body entropy distribution for equilib- 
rium state 

First, we describe the equilibrium state, which should 
be given by the state that maximizes the entropy under 
some constraints if the system is ergodic. 

Within the one-body approximation,^'^''' the entropy 
and the energy arc approximated by the one-body en- 
tropy Si and the mean field energy Hy^p using the coarse- 
grained density n{r) as. 



Si = 



j d^r n{r) log n{r). 



(18) 



Fmf = -\j d^r J d^r'n{r)G{r,r')n{r'), (19) 

where n(r) is normalized to N as in eq. (17). The angu- 
lar momentum / is evaluated by substituting the coarse- 
grained density n(r) into eq. (15). Then the theory be- 
comes simple and the maximum one-body entropy state 
is obtained® by maximizing eq. (18) with respect to n(r) 
under the constraints of N, Hmf a-nd / as 



n(r) = cxp(— 1 — a 



r — cr 



(20) 



with the Lagrange multipliers a, 6, and c. From eq. (8), 
satisfies 

V^(/)(r) = exp(-l - a + ^(r) - cr^), (21) 

which can be solved numerically . 

Note that Si would be a Casimir constant if the density 
field strictly followed the Euler equation (6). 



3. 2 Velocity distribution of randomly placed point vor- 
tex system 

The velocity distribution of particles is an important 
element in the kinetic theory. In an ordinary system, the 
velocities and the positions of particles are two sets of 
dynamical variables that define a state. In the present 
system of the point vortices, however, the velocities are 
determined by the positions, and its distributions has 
been analyzed in terms of Levy's stable distribution. '^"'"■^^ 

Let us consider the pancake density distribution, where 
the particles are located randomly with a constant aver- 
age density within a radius R. The system undergoes a 
rigid bulk rotation with the period r^ot on average. At 
the center of the system, the particles have no average 
velocity and only the velocity fluctuation exists. For large 
N, the velocity distribution at the center of the system 
has been found to be the Gaussian with the variance pro- 
portional to n log N for the small velocity region, and the 
power-law V^'^ in the large velocity region, 

' ;i™vexp(-;j^43iTv^')> iV«V,,,,{N)) 



W{V) = 



with the crossover velocity 

9 



(22) 



1/2 



V.rit{N)^ l^^logNj [log(logiV)]i/2, (23) 

where 7 is the circulation (or electric charge) of a par- 
ticle and normalized to unity in this paper. Note that 
the power law behavior of the distribution in large V 
comes from the distance distribution of the closest par- 
ticles among the randomly located particles. 

3.3 Kinetic theory of point vortex system 

Exploiting the similarity to the relaxation process of 
the self-gravitating system, Chavanis has developed the 
kinetic theory for the point vortex system. ^^'^^ As we 
have described, there are two relaxations, i.e. the fast 
relaxation and the slow relaxation; we will briefly review 
some of his results for the slow relaxation. 

After the system establishes a quasi-stationary state 
following the Euler equation, small effects of the collision 
term in cq. (11) gradually stir the system to cause the 
slow relaxation, or the coUisional relaxation. Chavanis 
has estimated the time scale of the coUisional relaxation 
in terms of the diffusion of point vortices caused by the 
collision term as follows. 

Assume that the velocity distribution is the Gaussian 
given in eq. (22) for the whole V region, then we can 
estimate the typical velocity V^yp from the mean square 
velocity as. 



"^typ 



(F^)^^ log TV. 



(24) 



For the pancake state with the flat distribution, there 
is no shear flow, thus the diffusion constant D may be 
estimated as 



D ^ /o^typ ~ 7 



(25) 



where Iq ^ l/y/no is the average inter-particle distance. 
For the state with singly peaked distribution, there 
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exists the shear flow around the peak, therefore, the dif- 
fusion is not isotropic. The diffusion coefficient for the 
radial direction may be estimated as 

7n(r) 



i?w is chosen as 



|S(r) 



■logA^ 



(26) 



where the correlation time r has been estimated as r 
1/|E| with the local shear 



|S(r) 



d fVir) 

' — 

dr \ r 



(27) 



If the slow relaxation is due to this diffusion process, 
the relaxation time Trciax would be the time for a particle 
to diffuse over the system size: 

L2 N 

Trelax ~ ~ — T==7"rot (28) 

for the pancake state without shear flow, and 



'^relax ^ , ^ 



N 



U log N 

for the singly peaked state with shear flow. In the both 
cases, Trelax Incrcascs almost linearly with N in units of 
the bulk rotation time. 

4. Simulation method 

We simulate the point vortex system by integrating 
eqs. (2) — (4).'^® Number of particles N ranges over 
N — 128 ~ 2048. Force acting on each particle is cal- 
culated by the simple sum of Coulomb force from all the 
other particles, therefore computational complexity of 
calculating the interactions is 0{N'^) for each time step. 
Numerical integration in time is performed using the 
second-order Rungc-Kutta method. The time step of the 
integration is At ~ 0.005Trot/47r. The error AH in the 
total energy H tends to increase almost linearly in time, 
and the average increase of the error is AH/H ~ 10~^ 
for TV = 128 and AH/H - 10"^ for N = 2048 per bulk 
rotation period Trot- 

Initial configurations of particles are generated using 
random number sequences to be consistent with a given 
macroscopic density distribution. The problem is that, 
if two particles are too close to each other, they rotate 
around each other with large velocity; This could cause 
large integration error. To avoid this problem, we restrict 
the distance between any pair of particles to be larger 
than a certain limit length proportional to ^jAt/N in 
an initial configuration. This restriction is implemented 
by re-generating a particle position if the newly placed 
particle is too close to the existing particles. It should 
be noted, however, that such initial condition does not 
ensure that the particles do not come close to each other 
in future. 

The unit length is defined so that the average square 
distance A'^ from the center to be 



2 _ X) ; ^ ? 



A" = 



N 



(30) 



which is the angular momentum (15) per particle, hence, 
is a constant of dynamics. The radius of the container 



2.9, 



(31) 



which is comparable to that in the experiments.®'^ With 
this value, the container docs not impose any conceivable 
restrictions on the relaxation process. 

5. Simulation results 

In this section we present our simulation results, 
mainly on the slow relaxation process in comparison with 
the kinetic theory by Chavanis. 

After describing a couple of system setups used in the 
numerical simulation (§5.1), we present general behavior 
of the system and show that the time scale of the slow re- 
laxation grows almost linearly with N (§5.2). The veloc- 
ity distributions are shown in §5.3. Particle trajectories 
are examined in §5.4 to find anomalous diffusion, which 
is analyzed in terms of the Levy flight by decomposing 
the trajectories into sequences of steps, whose length and 
duration distributions are found to be of the power laws 
(§5.5). 

5.1 Initial configurations 

We will present simulation data mainly for the system 
setups that starts from two different initial states, which 
we call Setup A and Setup B in the following. 

Setup A is the simulation that starts from the pancake 
state with a constant macroscopic density distribution in 
the region r < Rq, where the radius of the outer edge of 
the pancake Rq is determined by the condition (30). In 
the continuum limit, this pancake state is the minimum 
energy state for a given angular momentum, therefore, 
the system cannot relax macroscopically. 

Setup B is the simulation that starts from a double ring 
distribution. The particles are distributed in the regions 
0.4i?o < r < 0.6i?o and 0.8i?o < r < Rq. In this case, 
the system relaxes into a singly peaked state with shear 
flow, as we will see below. 

The numbers of particles N are same for both setups 
A and B. The relaxation properties of the macroscopic 
density are mainly obtained from Setup B of the double 
ring initial state. Setup A of the pancake initial state 
is used to investigate the diffusion process of individual 
point vortices. 

Another variation of initial configuration, which we 
call Setup B', is used in Fig. 6 in order to see the initial 
state dependence of the slow relaxation. Setup B' is the 
simulation that starts from a single ring distribution. The 
particles are distributed in the regions 0.4i?o < r < Rq. 
The energy happens to be almost the same with that of 
Setup B. The system relaxes into a singly peaked state. 

5.2 Fast and slow relaxation 

Figure 4 shows the time development of the one-body 
entropy Si{t) as a function of time t in the semilogarith- 
mic scale. The one-body entropy is calculated by coarse- 
graining the distribution of point vortices using a mesh 
shown in Fig. 3. The size of the cell at distance r from 
the center is dr = R^ /nir in the radial direction and is 
r56 = 27r i?w/™r va. the azimuthal direction, where is 
the number of division in the radial direction, thus the 
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shape of a cell is elongated in the azimuthal direction. 
This gives the area of a cell 6 A w {Sr){rS9) = 2 A/ ml 
where A = ttR^ is the total area of the system. Here, 
nir = 20 is used for all the simulations. 

We observe three stages in the development of Si{t): 
(i) diocotron instability {t < iTrot), where Si(t) is almost 
constant, (ii) fast relaxation (Iriot ^ t < 10T,.ot), where 
Si(t) shows a rapid increase, and (iii) slow relaxation 
(lOTi-ot ^ i)i where Si{t) increases slowly with time. 

The one-body entropy Si{t) as a function of time t 
in the fast relaxation is shown in Fig. 5 for various N. 
We do not observe the systematic A^-dependcncc of the 
relaxation time scale if we scale the time t by the bulk 
rotation period Trot oc 

Figure 6 shows the time evolutions of the density dis- 
tributions as a function of r. The initial states are the 
double ring initial state (Setup B) (a), and the single 




Fig. 3. An example of mesh used to observe the macroscopic den- 
sity field for one-body entropy calculation. Number of partition 
in azimuthal direction increases by one as r increases (in the fig- 
ure, 1,2 and 3 from the center to the periphery). In the actual 
observations, the resolution in radial direction m,r is = 20. 



ring initial state (Setup B') (b). The maximum one-body 
entropy state is also shown by the thick solid lines for 
comparison. The energy and the angular momentum are 
same for these two cases, so that the maximum one-body 
entropy state is same for the both cases. Although the 
quasi-stationary state achieved after the fast relaxation 
depends on the initial state, after the slow relaxation, the 
density distribution approaches the maximum one-body 
entropy state in both cases. Another example is shown in 
Fig. 7, where the initial single ring state relaxes into the 
vortex crystal state after the fast relaxation, but the vor- 
tices are smeared to merge into a single broad peak and 
the state eventually approaches the maximum one-body 
entropy state through the slow relaxation. 

Figure 8 shows that the time development of the one- 
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Fig. 6. Time development of the density distribution vs. r for the 
double ring initial state (Setup B) (a), and the single ring initial 
state (Setup B') (b). In the both cases, the number of particles 
is N = 1024 and the energy arc almost the same for the both 
cases. The maximum one-body entropy state is also shown by 
the thick solid curves. For each curves of simulation results, the 
density distribution is averaged over 100 samplings at different 
times within the duration of about 10 rotations. 




sim. A'=1024 - 
10 100 1000 

f/Tiol 



Fig. 4. Time development of the one-body entropy Si{t) in the 
semilogarithmic scale for Setup B with = 1024. 




Fig. 5. Time evolution of the one-body entropy Si during the fast 
relaxation for various N for Setup B. The time t is scaled by the 
bulk rotation time Trot oc 1/N. 
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Fig. 7. Time evolution of the point vortex distribution with = 
512. The initial state is the single ring state with a constant 
density at 0.9i?o < r < Rq of randomly placed particles. 
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Sim. ^=768 
Aexp(-f/(64T„,i))+S,, 



100 



Fig. 8. Time evolution of the one-body entropy Si during the 
slow relaxation for Setup B with N = 768. The simulation data 
are averaged over 100 successive points. The dashed curve shows 
an exponential decay of cq. (32), with the fitting parameters 



■■ 0.775, A = 0.19, and r^da 
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body entropy Si{t) during the slow relaxation can be 
approximated by the exponential form: 

Sit) w S-cq - Aexp(-t/r,.eiax). (32) 

We determine the parameters ^oq, A and Trdax by fitting 
with the simulation data. The relaxation time scale Tioiax 
is plotted against the number of particles N in Fig. 9. We 
find that Trdax increases almost linearly with N if Treiax 
is measured in the unit of the bulk rotation time r^ot cx; 
1/iV. This is consistent with the theory of Chavanis in 
§3.3. 

5.3 Velocity distribution 

Figure 10 shows the radial velocity distributions of the 
simulations of Setup A (the pancake state) (a) , and Setup 
B (the singly peaked state) (b). In both cases, we mea- 
sure the velocity distribution in the same time region 
SOTrot ^ t < 240Trot, which is within the slow relaxation 
for the case of the singly peaked state (b), while the pan- 
cake state (a) docs not show any macroscopic relaxation. 
Only the distributions of the radial component of veloc- 
ity are plotted because the azimuthal velocity contains 
the average bulk rotation. 

In Fig. 10(a), the Gaussian distribution that scales as 
v/ {N log N) is observed for small v, as is expected by the 
theory in §3.3.^^''^^"'^^ For large v, however, we find the 
exponential tail instead of the power law tail. 

The reason of the lack of the power law tail seems 
to be the restriction of the particle distance imposed on 



the initial state, as is explained in §4. If we remove this 
restriction and locate the particles randomly, then we ob- 
tain the tail closer to the power law as shown in Fig. 
11. Note that the precision of the energy conservation is 
poor for the simulation without the restriction. 

For the singly peaked state (Fig. 10(b)), we find that 
the distribution is broader than that for the pancake case 
(a) but the tail is still exponential. Although the distri- 
bution does not fit to the Gaussian very well, it is still 
scaled by v/{NlogN) in small v region. 

5.4 Anomalous diffusion 

Figure 12 shows a trajectory of the radial position of 
a particle in the quasi-stationary state. One can see that 
the motion of a particle looks like a random walk, while 
it occasionally undergoes long jumps with a step almost 
of the system size (Ar ~ 1). 

To see the diffusion behavior, the time development 
of the mean square radial deviation ((Ar(i))^) where 
Ar{t) = r(<)— r(0) is plotted in Fig. 13(a) for the pancake 
case with various TV (Setup A). We observe a superdiffu- 
sion 

((Ar(t))2) ^ t\ (33) 

with the diffusion exponent 7 = 1.75 ±0.1 > 1 for all A^. 
The square radial deviation saturates in ((Ar(i))^) w 0.1 
due to the finite system size L 1. From Fig. 13(b); we 
obtain 7 = 1.85 ± 0.1 for the case of Setup B (the singly 
peaked state). 

Figure 14 shows the "coefficient of anomalous diffu- 
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Fig. 9. Relaxation time Treiax of the one body entropy in the slow 
relaxation as a function of the number of particles A'. The j/-axis 
is normalized by the bulk rotation time r^ot. The simulation is 
by Setup B. 
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Fig. 10. Distributions of the radial component of velocity of two 
different simulation setups: Setup A (the pancake state) (a), 
and Setup B (singly peaked state) (b). The velocity is scaled 
by V*{N) = NlogN. The distribution functions are averaged 
over 200 samplings at different times during the slow relaxation, 

SOTrot < t < 240Trot . 
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Fig. 11. Distributions of the radial component of velocity for 
Setup A with (Q) and without (I) the restriction on the dis- 
tance between the particles (see §4). The number of particles is 
N = 128. The velocity is scaled by V*{N) = NlogN. The dis- 
tribution functions are averaged over 200 samplings at different 
times during the slow relaxation, SOrrot < < < 240rrot. 
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Fig. 12. Time development of the radial position of a particle for 
Setup B (the singly peaked state) with A'' = 1024. 
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sion" defined as 



Dn = 



(34) 



as a function of N . We find a power law dependence on 
N afiDN ^ N"^ with = 0.85 ± 0.1. 

Figure 15 shows the distribution function P(Ar|t) of 
the deviation Ar at various times t. The data are taken 
from the same simulation with that in Fig. 13 (a) for = 
2048. We see that the distribution curves overlap by the 
sealing, Ar/t^/"^ with 7 = 1.75. which shows the width of 
the distribution grows as V^'^. The tail of the distribution 
converges to zero faster than the exponential. 





Fig. 13. Time evolution of the mean square of the difference of 
radial position Ar(t) = r(t) — r(0) from its initial position in the 
log-log plot. The simulation setups are Setup A (the pancake 
state) (a), and Setup B (the singly peaked state) (b). 




Fig. 14. Coefficient of anomalous diffusion i?jv = {Ar'^)/t"' as a 
function of the number of particles A'^. The coefficient Djv arc 
estimated from the plots in Fig. 13. 7 = 1.75 is used for Setup A 
(the pancake state), and 7 = 1.85 is used for Setup B (the singly 
peaked state). 



f=0.080T„, 
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-1 -0.5 0.5 1 

Fig. 15. Distribution P(Ar\t) of the deviation of the radial posi- 
tion Ar{t) = r(t) — r(0) at various times t. The plotted data are 
same with those in Fig. 13(a) for N = 2048. The deviation Ar is 
scaled by with the same exponent 7 = 1.75 with Fig. 13(a). 



5. 5 Random walk and Levy flight 

In order to analyze the above diffusive trajectories as 
shown in Fig. 12, we decompose them into "steps". Each 
step is defined as the interval between two successive 
extrema in the radial motion Ar(t), thus each step has a 
step length I and a step time r. Then the radial position 
Ar and the elapsed time t at M-th step are expressed by 

M 



M 



Ar 



fc=i 



t 



(35) 



fc=i 



Figure 16(a) shows the distributions of the step length 
in the simulations of Setup A (the pancake state) with 
various TV, which exhibits the power law decay in the 
large step length. The step length distributions for vari- 
ous N seem to overlap with each other if we scale the step 
length as IVn, which indicates that the lower cutoff of 
the power law scales as Zmin ~ The upper cutoff 

of the power law decay is given by the system size, thus 
it does not depend on N. The power law distribution is 
expressed in the form of p(Z) ^ and the exponent 

/X = 0.7 ± 0.1 is obtained. We find that the distribu- 
tion function is almost symmetric both for the positive 
/ (jump towards the periphery) and the negative /(jump 
towards the center). 

Figure 16(b) shows the similar distribution function 
obtained for the step time t, which is normalized by the 
bulk rotation time Tiot oc 1/A^. The power law exponent 
X is defined by p{t) ^ t^'-^^^^ and x = 0.55 ± 0.1 is 
obtained. The lower cutoff of the power law scales as 

Tmin ~ 1/^. 

Figure 17 shows the step length distribution p{l) and 
the step time distribution p{t) for the singly peaked state 
(Setup B). 

Table I summarizes the exponents obtained in our sim- 
ulations; fi^s and x's are the power law exponents of the 
step size distributions p{l) and p(t), respectively, and 
their values are obtained from Figs. 16 and 17. The ex- 
ponents of the anomalous diffusion 7 from the simulation 
results in Fig. 13, and the exponents 1] for the N depen- 
dence of the anomalous diffusion coefficient Dn ^ iV 
from Fig. 14 are also listed. The exponents 7cst and 7;cst 
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Fig. 16. Distribution functions for the step length I for the radial 
direction scaled by \/~N (a), and step time r* = r/rrot in units of 
the bulk rotation time Trot (b). The simulation is by Setup A (the 
pancake state). The step length I is defined as a distance between 
two neighboring turns in the trajectories of the radial position 
of a particle, and t is duration time needed for a step. For each 
simulation with TV, 64 particles are randomly chosen and their 
trajectories arc analyzed in the slow relaxation, SOrrot < < < 
240Trot. 
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arc the exponents derived from the step distribution ex- 
ponents (see §6). 

Now we examine the correlation in the step sequences. 
Let the step length and step time of the fc-tli step be Ik 
and Tfc, respectively. The correlations {loh) and (roTfe) 
are plotted against k for the pancake state (Setup A) in 
Fig. 18. These correlations show the exponential decay, 
which indicates that the correlation in the step sequences 
is short-range. 

In order to analyze the step sequence in more detail, we 
define P(|Ar| | Af ) and P{t\M) as the distributions of the 
distance Ar and the time t at the Af-th step (eq. (35)), 
respectively. Contour plots of these distribution functions 
are shown in Fig. 19. The value of the distributions are 
shown in the proportion to the peak value Ppcak(A/) of 



(b) 




0.01 1 100 

radial step length Nn 




flight time t^=z/ Tro, 



Fig. 17. Distribution functions for the step length I for the radial 
direction scaled by \/7V (a), and the step time t* = r/riot in 
the unit of the bulk rotation time Trot (b). The simulation is by 
Setup B (the singly peaked state). The data are collected in the 
same way as in Fig. 16 



Table I. Exponents for the stop time distributions x, the step size 
distributions fi, and the anomalous diffusion exponents 7 and rj 
and the derived exponents 7cst = 2x7^ rjost = 2x7/^ — 1 from 
X and fj, See texts for the definition of each exponent. 



Setup 


X 




7 


V 


A 


0.55 ±0.1 


0.7± 0.1 


1.75 ±0.1 


0.85 ±0.1 


B 


0.55 ±0.1 


0.6 ± 0.1 


1.85 ±0.1 


1.2 ±0.1 


Setup 






7ost 


Vest 


A 






1.6 


0.6 


B 






1.8 


0.8 



P(|Ar| |Af) or P{t\M) for each step Af so that one can 
compare the width of the distributions between different 
M's. We observe that the width of P{\Ar\ \M) and the 
peak position of P{t\M) increase by the powers of A/, as 

(Ar)2 ^ M^/^'\ t - M^/^" (36) 

with n* = 1.1 ± 0.2 and x* = 0.85 ± 0.1. 

Note that the contour curves are almost equally spaced 
in the logarithmic scale, which means that the distribu- 
tion functions have the power-law tail. The growth of the 
tail of the distribution function P(|Ar| \ M) is saturated 
for large A/ where the tail approaches Ar ~ 1, or the 
system size. 

6. Analyses 

As we have shown in §5.5, the step length distribution 
p{l) and the step time distribution p{t) have the power 
law tail, and the correlation between the steps decays ex- 
ponentially. Based on these observations, we analyze the 
diffusion process in terms of the Levy flight and derive 
the exponents. Some of the basic properties of the Levy 
flight are summarized in Appendix. 

The typical radial deviation of position Ar^f and the 
typical elapsed time Im after the A/-th step can be esti- 
mated by Levi flight as 

i\± 'mini 



Ar^ 



, 72 



(37) 



from eq. (A-8). 

These should be compared with eq. (36), which are di- 
rect observation of P{Ar\M) and P{t\M). Both of them 
show the power law of Af , but the values of the expo- 
nents n* and X* estimated from the direct observation 
are not very close to the values of n and x from the step 
distribution, we suspect this discrepancy is due to the 
finite size effect as we have mentioned at the end of §5.5. 

Keeping this limitation in mind, let us continue the 
analysis based on the Levy flight picture; The expression 
(37) show that the total sums of eq. (35) are of the same 
order of the largest step out of the Af steps. This means 
that most of the contribution to eq. (35) comes from the 
largest step. We assume this largest step to be the same 
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Fig. 18. Correlation functions in the step sequences of the step 
length {lolk) (^)i the step time (toTj.) (b) as a function of k. 
The simulation is by Setup A (the pancake state) with N = 2048. 
These data are obtained from trajectories of randomly chosen 64 
particles in the slow relaxation process, 80rrot < i < 240Trot. 
The fitting parameters are C; = 0.00107 and Ct = 0.539. The 
inset in (a) shows a linear plot of stop correlation in the same 
range of k. 
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Fig. 19. Contour plots for the distribution of radial position 
P{\Ar\ \M) (a), and that of the time P(t*\M) (b) as a function of 
the number of stops M in the logarithmic scale, t* = t/Trot and 
the distribution functions are normalized by the peak distribu- 
tion Ppeak(A^) for each step M. The thick lines are | Ar| oc M^/ 
(a) and t* oc M^/^" (b) with fi* = 1.1 and x* = 0.85 that 
are obtained by fitting. The data are obtained from trajec- 
tories of randomly chosen 64 particles in the time duration, 
80rrot < t < 240rrot . For both figures, the simulation is by Setup 
A (the pancake state) with Af = 2048. 
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step for A?' and t. This is natural assumption because 
shorter steps are results of waggling motion caused by 
nearly vortices, vortices tend to move faster during the 
shorter steps. 

Using the assumption, the typical radial position Art 
at time t can be estimated by simply eliminating M in 
eq. (37). Identifying as t, we have 



Ar? 



with 



D 



N 



?2 . 

mm 

2x/m ' 



(38) 



(39) 



Comparing eq. (38) with eq. (33), we see that the anoma- 
lous diffusion exponent 7 is related with x and fi as 



7 = 2x/Ai- 



(40) 



The dependence of the anomalous diffusion coefficient 
Dn on the number of particles N is also obtained from 
eq. (39). Assuming that the short length cutoff of the 
step length distribution /min is approximately given by 
the average particle distance Iq, thus we have 



It 



N 



(41) 



Similarly, we assume that the short time cutoff of the 
step time distribution Tmin is given by the rotation time 
of the vortex pair tq with the distance Iq: 



To 



(42) 



where vq is the rotation velocity of the paired vortices 
with the distance Zq: and given by 



q 

vo ^ — , 
'0 



(43) 



from the 2-d Coulomb law. Here, q is the charge of a 
particle, which is normalized to unity. Note that these 
N dependences are consistent with the scaling form of 
the step size distributions p{l) and p{t) in Figs. 16 and 
17. Substituting these cutoffs into eq. (39), we obtain N 
dependence of the coefficient Dpf, 



D 



N 



thus the exponent 77 as 

V = 2x//^ -1 = 7-1. 



(44) 



(45) 



From the simulation results and the definition of the 
exponents, we obtain x = 0-55 and fi = 0.7 for the radial 
direction in the pancake state (Setup A). Thus, from eqs. 
(40) and (44), one obtains 7 « 1.6 and 77 « 0.6. Simi- 
larly, we obtain the exponents for the singly peaked state 
(Setup B). They are hsted in Table I. 

7. Summary and discussion 

We have performed the numerical simulations on the 
two-dimensional point vortex model with a unit circula- 
tion of the same sign, in order to study the relaxation 
process of pure electron plasma under the strong mag- 
netic field. Due to the long-range interaction between the 
vortices, the system behaves very differently from ordi- 
nary systems with short-range interaction. 



We have found the following: (i) There exist two types 
of relaxation: the fast relaxation and the slow relaxation. 
The fast relaxation takes place with the time scale com- 
parable with the bulk rotation time Trot and leads the 
system to a quasi-stationary state following the Euler 
equation, while the slow relaxation takes place after the 
fast relaxation with the relaxation time Treiax ~ Nr^-ot 
due to the individual motion of vortices; This is consis- 
tent with Chavanis' theory except for the logiV correc- 
tion. After the slow relaxation, the system reaches the 
maximum one-body entropy state, (ii) Individual mo- 
tion of point vortices in the slow relaxation process is 
superdiffusive with the exponent 7 w 1.75 for the pan- 
cake (shear free) state and 7 w 1.85 for the singly peaked 
(with shear) state. The coefficient of anomalous diffusion 
depends on the number of particles N in the power laws, 
(iii) The superdiffusive motion of individual vortices can 
be decomposed into a sequence of steps. The correlations 
of the step length and that of the step time along the se- 
quence are short range, and the distributions of the step 
length and the step time are of the power laws with the 
exponents fi and respectively. The superdiffusive mo- 
tion can be reconstructed from the Levy flight, i.e. the 
exponent 7 for the anomalous diffusion and the exponent 
77 for the N dependence of its coefficient are expressed 
as 7 = 2x//i and 77 = 2%/^ — 1; the former expression 
agrees with the exponent from the simulation within the 
error bars, but the agreement of the latter is not good. 

Among these results, our result (i) of the N depen- 
dence of the slow relaxation time seems to agree with 
Chavanis' estimation. However, his picture that the re- 
laxation is due to the normal diffusion of the point vor- 
tices is not confirmed by our simulation results, but we 
observe the superdiffusive behavior in the simulations. 



If we estimate the relaxation time r* 



with anomalous 



diffusion by the similar way as he did for normal diffu- 
sion, then 



clax 



Dn 



1/7 



N 



l-i)/7 



■''rot ''"relax ^ 



(46) 

which gives much shorter relaxation time than that ob- 
served. This implies that the motions of point vortices 
are not independent of each other and provide only weak 
mixing. Actually we observe in our simulations that a 
point vortex tends to avoid "collisions" during the long 
jumps. 

In the literature, ^^'^^ it has been pointed out that a 
test particle is convected for a long way comparable to 
the system size due to the existence of long-living large 
vortices. The diffusion coefficient proposed by Taylor and 
McNamara depends on time and converges to a constant 
only when (Ar/L)^ ^ 1 holds, which suggests anoma- 
lous diffusion. Although these results are for the neutral 
plasma, its behavior of the anomalous diffusion seems 
consistent qualitatively with our simulation results on 
the non-neutral plasma. 

Recently, Dubin and Jin performed fairly large-scale 
simulations on the 2-d point vortex system with a posi- 
tive charge, and determined the diffusion coefficients 
in the states without mean shear. Their results of the 
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diffusion coefficient show the A^^^^ dependence, which 
is consistent with that expected by Taylor and McNa- 
mara in the converging hmit, but their vahies of diffusion 
constant from the simulations seem substantially smaller 
than those predicted by the theory: This may be due to 
the slow convergence mentioned in the above paragraph. 
They have also shown that the diffusion in the radial 
direction is normal under the external shear. 

In experiments, the number of electrons is very large 
while the charge of each electron is very small. In the 
present model, the situation may correspond to the case 
in the limit of the infinite N with a fixed total charge, 
namely, a fixed bulk rotation time Trot, in which limit 
the slow relaxation never takes place and we only observe 
the quasi-stationary states. If this is the case, actual slow 
relaxation that may be observed in experiments should 
be due to non-ideal effects such as three dimensionality of 
the Malmberg trap, scattering at the end of the cylinder, 
impurities, etc. 
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Appendix: Levy flight and anomalous diffusion 

In this appendix, we summarize some basic formulas 
of the Levy flight. 

Consider a random walk problem where the size of each 
step I is a random variable without correlation, then after 
the M-th step the particle position X is 



distribution: 



M 



(A-l) 



If the distribution p{l) of step size / has a finite second 
moment (Z^), the central limit theorem tells us that the 
distribution of X at the A/-th step is known to be the 
Gaussian distribution for large M and the dispersion in- 
creases as 



{{X-{X)f 



DM, 



(A-2) 



where D is the diffusion coefficient. This represents the 
normal diffusion law. 

If the distribution function p{l) is the power law with 
the diverging second moment. 



(A-3) 



for large |/|, then the asymptotic form of the distribution 

of X at large M is given by the Levy's stable distribu- 
tion.41.42 



1 



P{X\M) 



1 .L ( ^~^^y 



(0<M< 1), 
(1<M<2), 

/A-4) 

where L^^p{x) is the scaling form of the distribution of 
X &i M ~> oo. The parameter (3 characterizes the degree 
of asymmetry, and is determined from the asymmetry of 
p{l) for large \l\. The /? = case represents the symmetric 



L^,q{Z) = ^ r dfcexp(zfcZ-qA:r), 



(A-5) 



where C is a scaling factor. This gives the Gaussian dis- 
tribution when ^ = 2. As for the /3 = 1 case, the distri- 
bution is given by 

dsexp(sZ- C's^), (A-6) 



d—io 



which is zero for Z < representing a completely asym- 
metric case. The distribution L^^p(Z) has an approxi- 
mate form of L^^p{Z) ~ ^-(1+'') for large Z. 

The fact that the above distribution function P{X\M) 
scales as XjM^I^ can be understood"'^ by introducing 
the effective cutoff Zc, to the step length distribution 
By simple argument, one can see the cutoff depends on 
the number of steps M as 



/e - Af^/^U„, (0 < M < 1) 
Then, a typical walk may be estimated as 



(A-7) 



M 



Xm = Y.^^~ 



dllp{l)^ M'/f^Un, (0<M<1), 

(A-8) 



which means the distribution function scales as X/M^^f^. 
Therefore, peak position Xp ~ M^^^^ and the square of 
typical width X^ ~ Af^/^ of the distribution P{X\M) 
at the Af-th step can be considered as measures of the 
anomalous diffusion by the Levy flight, even though the 
distribution function has diverging moments. 
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